ADAPTIVE NUMERICAL INTEGRATION AND CONTROL VARIATES FOR 

PRICING BASKET OPTIONS 
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Abstract. We develop a numerical method for pricing multidimensional vanilla options in the Black-Scholes 
framework. In low dimensions, we improve an adaptive integration algorithm proposed by two of the authors by 
introducing a new splitting strategy based on a geometrical criterion. In higher dimensions, this new algorithm is 
used as a control variate after a dimension reduction based on principal component analysis. Numerical tests are 
performed on the pricing of basket, put on minimum and digital options in dimensions up to ten. 
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1. Introduction. Pricing financial derivatives generally boils down to numerically compute 
an expectation, at least for European type contracts. Basically, there are two ways of doing so: 
either solving a partial differential equation or resorting to Monte Carlo techniques. Monte Carlo 
integration is known to provide better results when the dimension of the problem increases but 
for reasonably small dimensional problems, its efficiency is not that clear. However, alternative 
numerical integration techniques may help in such situations. In the Black-Scholes framework, 
the pricing of vanilla options reduces to a numerical integration problem over Mf^, where d is the 
number of underlying assets. This problem has some specificities coming from the properties of 
the function to integrate. First, the integrand decreases quickly away from the origin and hence we 
can consider that the integration domain is a hypercube [—A, A]"^ with A relatively small. Second, 
the integrand is clearly not a smooth function in the whole domain as it is only continuous at 
the interface between the area where the function vanishes and the area where it is positive. 
Moreover, this interface is neither known and nor located at the boundary of the hypercube so 
there is no hope of using techniques like the periodisation method [7j to increase the smoothness 
of the integrand. Finally, the integrand is a muldimensional function in a dimension d, which can 
be large, so one also has to deal with the curse of dimensionality. 

Hence, it seems quite natural to use Monte Carlo or quasi-Monte Carlo ]^7\ methods to face 
such a difficult problem. In fact a crude use of these methods is not necessarily sufficient to reach 
a good accuracy. Besides the usual variance reduction methods, one needs to develop adaptive 
methods to make them really competitive. In many situations, the function to integrate or to 
approximate may have completely different behaviors in terms of variations or even in terms 
of regularity in different parts of the domain D. In those cases, it might be more efficient to 
adaptively split D in subregions according to error indicators based on quadrature points. The 
most famous adaptive Monte Carlo integration routines are MISER [5D] and VEGAS [12] ■ They 
rely on stratified and importance sampling and error indicators based on the empirical variance, 
respectively. Quasi-Monte Carlo versions of these two algorithms have been introduced in [21]. 
More recently, adaptive approaches have been developed for stratified sampling [SI |S] and for 
importance sampling [11] . The idea of adaptive algorithms is to make use of past simulations to 
help better leading future ones. In the case of stratified sampling, both the number of samples in 
each strata and the boundaries of the strata can be learnt inline. The idea remains quite similar 
for adaptive importance sampling, which consists in adaptively learning the optimal change of 
measure from the already drawn samples. Following the methodology of [22], two of the authors 
have developed an adaptive integration algorithm [3] based on quasi-Monte Carlo quadratures, 
which proved to be efficient for very smooth but also for less smooth functions [TB] • This algorithm 
has already been tested on the pricing of basket options achieving excellent results in dimension 
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two but loosing most of its efficiency for larger dimensions. The loss of accuracy is mainly due to 
the difficulty to capture the interface between the two regions of interest. 

Our goal is to improve this algorithm in order to use it in higher dimensions, for the pricing 
of more general options like digital or put on minimum options and also for the computations of 
sensitivities with respect to the parameters of the Black-Scholes model. The improvements of the 
algorithm will consist in a new splitting criterion and on a dimension reduction using a principal 
components analysis combined with control variates. 

The rest of the paper is organized as follows. In section 2, we describe the different kinds of 
options considered in the Black-Scholes framework and we discuss the type of sensitivities we are 
interested in. In section 3, we recall the adaptive algorithm developed in [3\ and introduce a new 
criterion based on very simple geometric considerations like the ones used in mesh refinement for 
finite element methods. In section 4, we compare the new criterion to the old one on the pricing 
of different types of options in low dimensions. We focus on examples in dimension two in order 
to emphasize the quality of the mesh refinement. Section 5 is devoted to sensitivity analysis. 
We compute the Delta of the different derivatives using standard deterministic techniques like 
polynomial interpolation allowed by the high accuracy of our adaptive pricing method. The last 
section deals with higher dimensional models. The adaptive method is coupled with dimension 
reduction via a principal component analysis to derive a variance reduction method based on 
control variates. 

2. Presentation of the model. 

2.1. The model framework. We consider a Black-Scholes model in dimension d in which 
each asset is supposed to follow the standard one dimensional dynamics given under the risk 
neutral measure by 

dSl = Slirdt + aidWi) 

with S'q = and where Wt = (VF*, M^j) denotes a vector of correlated standard Brownian 
motions. The volatility cr is a vector in Mf^, the instantaneous interest rate is r and {s^,...,s'^) 
is the vector of spot values. The covariance structure of these correlated Brownian motions is 
supposed to be defined by {W, W)^ = Tt where F is a positive definite matrix with all its diagonal 
terms equal to one. In the numerical examples considered in the next sections, we assume that 

Tij = 5i^j + p{l - S,j) 

where the parameter p e] — jzjtM to ensure that the matrix T remains positive definite. We 
introduce the Cholesky decomposition C of F (such that CC^ — F) and denote by Ci its z*'* row 
for 1 < i < d. The Black-Scholes model can be rewritten 

where now i? is a standard d— dimensional Brownian motion. 

In this model, we want to price options with payoffs written as functions of the asset price at 
a maturity time T. Hence, the price is given by the discounted expectation exp(— rT)E(?/;(5'T)) 
where the function tp characterizes the option type. In the following, we consider three different 
multidimensional options for 

• Basket options with payoffs 



• Digital basket options with payoffs 
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• Put on minimum options with payoffs 



V'(5t) ^ [K- min 5* 



) 



+ 



The vector (Ai, Xd) represents the weight of the different assets within the basket. These weights 
may be negative to allow to consider exchange options. The variable K denotes the strike price 
and the vector U corresponds to an upper barrier on the asset price at maturity time T. Most of 
the time, the expectation E{tp{ST)) needs to be computed numerically by means of Monte Carlo 
methods. In order to use a systematic approach, the expectation E(i/)(5t)) is usually rewritten 
as E((/)(G')) where G is a standard normal random vector in M.'^ and (/) : K** — >■ K"*" is a measurable 
and integrable function. To use our adaptive method described in the next section, we need to 
transform the computation of E(0(G)) into the numerical computation of 



where p{x) is the density of G and A G IR+ is chosen large enough to have a good approximation 



2.2. Key role of the delta. When selling a financial derivative, the price of the option 
is obviously of a primary interest but one should keep in mind that the original definition of an 
option price goes back to the replicating theory. The price is actually defined as the value at time 
t = of the replicating portfolio. Hence, the price becomes fairly useless if we do not know how to 
implement the replicating portfolio; hopefully, we precisely know how many stocks the portfolio 
should carry on at any time t and this quantity is given by the famous Delta of an option defined 
as the gradient of the option price with respect to the spot vector 



There are several numerical methods for computing the delta of an option. The most commonly 
used method is based on a finite difference approach because of its automatic application. However, 
when coupled with a Monte Carlo method it often yields a poor accuracy unless a very large number 
of samples is used. 

3. Description of the adaptive algorithm. 

3.1. Quasi-Monte Carlo quadratures. Any adaptive integration method relies on a 
quadrature rule designed for the non-adaptive case. While Monte Carlo or quasi-Monte Carlo 
methods can deal with integrands with no or little regularity, usual quadrature rules like Gauss 
product rules are built for functions having a given regularity or belonging to particular bases. 
These quadratures should not be too sensitive to the dimensional effect and thus quadrature for- 
mulae based on interpolation or approximation on different reduced size bases are considered. In 
the case of c?— dimensional Fourier bases on periodic smooth functions, Korobov spaces [S] are 
built using that the Fourier coefficients am verify 



where fh = max(l, |m|) and /3 > 1 is linked to the regularity of the integrand. The corresponding 
quadrature formulae are lattice rules O (2^ which are exact for the most significant Fourier 
coefficients according to this decay. Non-periodic but smooth functions can be periodized T to 
still use these quadratures but with an increasing constant of decay. 

In the case of polynomial approximations, Novak and Ritter |18j have obtained quadrature 
formulae, which are exact for polynomials with total degree less than a given value. Based on the 
use of the control variates method on piecewise interpolation polynomials, Atanassov and Dimov 
[T] built a numerical method reaching an optimal rate of convergence for multivariate smooth 
functions with a fixed degree of differentiation. 




ofE(0(G)) by /(A). 



A^\/soexp{-rT)E{yj{ST)). 
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The previous quadratures give highly accurate results for smooth functions but they can yield 
a really poor accuracy for non-smooth functions. Monte Carlo methods are somehow more stable 
because they are not sensitive to the smoothness of the integrand. The quadratures developed 
in [15] combine the approximation on reduced Tchebychef polynomial basis and the use of quasi- 
Monte Carlo points to build this approximation. They are especially efficient for very smooth 
functions but they can also handle pretty well only continuous functions. They have been obtained 
after successive improvements of an initial adaptive Monte Carlo method |13j via quasi-random 
sequences [12] and the introduction of Tchebychef polynomial basis of Korobov type [13]. We 
recall these formulae for a multivariate function defined on the hypcrcube [—1,1]''. For m S N, 
let TO = max(l,TO). We introduce the set 

Wd,q = {m e {mi...md) < q} 

which corresponds to the level q of approximation. The reduced Tchebychef polynomial approxi- 
mation writes 

meWd,q 

where the L^.q — card{Wd,q) coefficients bm verify 



\b„A < 



{mim2...md 



\2L 



for a C^^ function. To compute the numerical approximation of the coefficients bm, we fit the 
model 

bmTmiixi)Tm2{x2)-Tma{xd) 

meWd,q 

to the observations of the function / at some quadrature points Pi = {X^^\ X^'^'') with 1 < * < 
M. The choice of the points Pi is obviously crucial for the condition number k{A) of the least- 
square matrix A. We have proved in [16] that if these points are independent random variables 
distributed according to the multidimensional Tchebychef density 

w{xi,X2, ■■,Xd) = Yi — / 

i=i '^V 1 ~ 

then k{A) goes to 1 when M — >■ cx) at a Monte Carlo speed Moreover the variance is 

bounded by one for any set Wd^q- An even better choice to represent the density w is to use 
quasi-Monte Carlo points or quantization points [19j . 

We compute numerically the inverse of the matrix A in order to obtain once and for all 
quadrature formulae for each of the coefficients bm, and finally a quadrature formula for the integral 
itself. In the following, the quadrature points are built using a x Ld^q Halton points with one 
additional point at each corner of the domain for a total of M = a x Ld^q -\- 2"^ points. In most 
situations, the parameter a is chosen equal to 3 which is sufficient to ensure a small value for k{A). 
The corner points are control points to detect a possible change of regularity of the function /. 
Nevertheless, we will see in Section 4, that more quadrature points may be necessary to detect 
this change of regularity on the difficult example of digital options. The approximation of 

/(/) = / J{x)dx 

is given by the quadrature formula 
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which is obtained thanks to the corresponding approximations Qd.a,q,b„Sf) of the coefficients 6„i 
and the integration of the approximation model. We also denote by Qd,a,q,B{f) and Qd.a,q,b^,Rif) 
the relative quadrature formulae on a given rectangle R which will be used in the adaptive inte- 
gration algorithm. 

3.2. Error indicators. We keep the same error indicators than the one used in our previous 
paper 3 . They rely on hierarchical quadratures for the integral of / but also for some of the 
coefficients of its weighted mean-square approximation. We select two sets Wd,qi and Wd,q2 with 
1 < qi < 92, the corresponding quadrature formulae Qd.,a,qAf) and Qd,a.,q2{f) for the integral of / 
and also the quadrature formulae for some of the leading coefficients in the approximation model. 
These coefficients are the d+1 coefficients belonging to the set Ad of the coefficients 6m for which 
all the indexes (toi, TO2.., rrid) are equal to zero or only one is non-zero and equal to one. Our error 
indicator £^d,a,gi. 92 ,-«(/) is 

\Qd ,a,qi ,R if)- (/)l 

on a given hyperrectangle R. This indicator is more robust than the usual indicator based on the 
comparison between only Qd.a,qi,R{f) and Qd,a,q2,R{f)- Indeed, it is less likely to happen that all 
the estimators of the leading coefficients are close to each other but all wrong. The approximate 
value of the integral on a given hyperrectangle is Qd.a,q2,R{f) and the approximate integral of / 
is the sum of the integrals over all the hyperrectangles. 

3.3. Splitting strategies. Now, we describe the splitting strategies used in the numerical 
experiments. At each step of the algorithm, the list of all the remaining hyper-rectangles involved 
in the algorithm is stored and these hyperrectangles are sorted according to their error indicator. 
Once the splitting is performed, the hyperrectangle R with the largest error indicator is removed 
from the list and the hyperrectangles Ri and R2 are inserted in the list according to their error 
indicators Ed,a.,qi,q2.,Ri{f) and Ed,a,qi,q2,R2 (/)■ 

3.3.1. Fully adaptive splitting. The first strategy named FAS is fully adaptive. It consists 
in trying the d possible ways to divide the hyper-rectangle R with largest error indicator in two 
equal size pieces Ri and R2 along one of the axes and keeping only the best splitting. The best 
splitting is the one for which 

Ed ■a,qi -m-.Rl if) + Ed ,a,qi ,g2,-R.2 (/) 

is minimum among the d possible choices. Each step of this algorithm requires dN{aLd.q -f 2'^) 
evaluations of the function /. The computational cost of the method should also take into account 
the cost of inserting the two newly created hyperrectangles into the list of all hyperrectangles. 
This can be done efficiently by an insertion sort which has a complexity of 0{N \og{N)) for our 
algorithm with N steps. In most practicals examples, the evaluation of the function / requires a 
lot of atomic computations and hence it makes sense to neglect the cost the sort compared to the 
number of function evaluations. 

3.3.2. Geometrical Random Splitting. The FAS requires d trials to find the optimal 
splitting. It is mainly interesting when the integrand has variations different from several degrees 
of magnitude from one coordinate to another like for instance for the two dimensional function 
cos(200a; -I- y). When pricing vanilla options, this kind of situation is unlikely to happen. Further- 
more and maybe more importantly, we have noticed in our previous work [3] that the convergence 
problems of our algorithm came from a too fine splitting in one direction near the interface between 
the regions where the regularity changes. 

We propose another strategy named geometrical random splitting (GRS) in order to reduce the 
computational costs and to solve the convergence problems by finding more precisely this interface. 
At each step of the algorithm, we divide the hyperrectangle with the largest error criterion in two 
equal size pieces i?i and R2 uniformly at random among the admissible directions. The admissible 
directions are the ones having the larger length among the axis of the hyperrectangle. The initial 
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hyperrectangle is always a hypercube which means that we can never have a ratio more than two 
between the different lengths of the axis of the hyperrectangles occurring in the mesh. The cost 
of the algorithm is N{aLd.q + 2^^) and it is stochastic which enables to compute some statistics on 
its results. 

4. Pricing Vanilla options in low dimensions. In this section, we compare the FAS and 
the GRS for pricing vanilla options on several examples already treated in [3] and also on examples 
from the general model. On all numerical examples, the same quadrature formulae will be used 
for the two methods. The number of iterations of the algorithm will be = 2000 for the FAS 
and consequently d times more for the GRS in order to keep the same complexity. 

4.1. Basket options. Let us focus a little on basket options. The price of such options, with 
payoffs only depending on the asset price at maturity time T, can be expressed as an expectation 
e~''^E(^/'(S'T)). When the random vector St has a density fs with respect to the Lebesgue 
measure, computing this expectation is easily turned into a numerical integration problem 



E(^(5t)) = / iP{x)fs{x)dx. 

J [— C30,00]'' 

When = 1, there is no need of a numerical integration as a closed formula exists. In the 
Black-Scholes framework, we know that 

St (^S'o exp | (r - + a.QVf 

where G is a random normal vector with values in M.'^. Hence, 

Eii^iST)) = ^ (^S'o exp |(r - + a.QVT x'^ ,l<i<d^ _L_e-l-lV2rf^. 

For the particular case of a call basket option, the payoff ip writes down 

V'(s) - (j2 - ^) ■ 

Therefore, if we denote by V{T, K) the price of the call basket option, we can write 

Using the equality (s — K)^ — {K — s)_|_ = s — K and letting the price of the corresponding put 
basket option be U{T, K) — e^''-^E ^(AT — X^iLi ^iST)+j ; we obtain the so-called call put parity 
relationship 

d 

V{T,K) - U{T,K) = - A:exp(-rr). 

i=l 

In general, this formula is used to obtain the hardest price to compute (in terms of variance) 
between the call or the put option prices from the other. Here, we independently compute the 
call and put option prices and use this formula as a criterion of accuracy. The infinite integral is 
truncated and we let 

V{T,K,A) = e-^'^j^^j^ ^^^^ |^^A,5Sexp|(r-^)r + a.QVrx|-xj e-N^/^dx 

be the truncated estimation of V{T,K). The truncated approximation of the put basket option 
price U{T, K, A) is defined in a similar way. Note that these multi-dimensional integrals are 
truncated on a square domain. In all the following examples dealing with basket options, we fix 
the weights of the basket as Ai = ^ so that the baskets are all homogeneous and their weights sum 
up to one. 
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4.1.1. Examples in dimension 2. The numerical tests considered in this paragraph are 
very similar to the ones treated in [5]. Three different values for the strike price are tested, one 
at the money Ki = 100, one out K2 = 127.80 and one completely out = 300. These examples 
were already studied in and the conclusion was that the FAS was outperforming Monte Carlo 
and quasi-Monte Carlo integration. Hence, we only compare the GRS to the FAS. Moreover, the 
new algorithm is only run once as we have a deterministic error criterion based on the call-put 
parity formula. 

The GRS is used with quadrature formulae of degrees 18 and 24. The latter corresponds to a 
number of function evaluations of 2000 x 403 x 2 ~ 1.6 x 10^. We give in Tables HTTI and W7^ the 
values of the truncated estimations for Ai = 12 and A2 — 13. We also compute 



C(T, K, A) 



1 



V{T, K, A) - U{T, K, A) - -S^ 



,(1) 



1 



for these reference values and denote by Coid{T, K, A) the same error indicator using the FAS. 
Truncating the domain of a real valued standard normal random variable to [—12, 12] may seem 
far too large. This is true from a Monte-Carlo point of view and we have actually also run some 
tests with more conventional values such as A = 5. These tests were already showing a far better 
accuracy than a crude Monte Carlo; the results obtained on the examples of Table 14.11 already 
had 4 accurate digits. However, since we are targeting to compute sensibilities, we need a better 
accuracy, which explains the choices of the parameter A. 





V 


U 


C 


Cold 


{Ki,A,) 


28.49407706 


14.564874729 


2 X 10-** 


1 X 10-« 


{Ki,A2) 


28.49407708 


14.564874726 


2 X 10-1" 


1 X 10-** 


{K2,A,) 


18.85549194 


28.853971355 


2 X 10-** 


2 X 10-** 


{K2.A2) 


18.85549196 


28.853971353 


2 X 10-** 


7 X lO-'J 


(i^3,^l) 


1.810536572 


160.02292952 


2 X 10-** 


2 X 10-** 


{K^,A2) 


1.810536593 


160.02292952 


4 X 10-1" 


4 X 10-1" 



Table 4.1: Basket option with parameters T = 3, 
5^^^ = 5^^^ = 50, (Ti = 0-2 = 0.4, p = 0.3 



0.05, d 



On this first example, we obtain a very good accuracy of at least 8 digits on all the computa- 
tions of the prices of call and put options. These values are both validated by the parity call put 
formula and by the comparison between the truncated approximations for the two different values 
of A. There is no significant difference between the two splitting strategies. The two versions of 





V 


U 


C 


Cold 




20.04091112 


6.1117087694 


2 X lO-'J 


1 X 10-** 


{KuA2) 


20.04091112 


6.1117087676 


1 X 10-1" 


1 X 10-** 


{K2,Ai) 


8.915343209 


18.913822596 


7 X 10-1" 


2 X 10-** 


{K2.A2) 


8.915343211 


18.913822598 


5 X 10-1" 


7 X lO-'J 


(i^3,^l) 


0.021755879 


158.23414880 


4 X 10-1" 


2 X 10-** 




0.021755880 


158.23414880 


6 X 10-11 


4 X 10-1" 



Table 4.2: Basket option with parameters T = 3, r = 0.05, d — 2 
5'(i) = S'^Q^ ^ 50, 0-1 = 0-2 0.2, p = 0.7 



the algorithm are still very efficient on an example with a smaller volatility and more correlated 
assets. The accuracy is even better than on the first example and we note that the GRS is now one 
or two digits more accurate than the FAS. It is also interesting to compare the meshes obtained 
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(a) Geometrical Random Splitting (b) Fully Adaptive Splitting 

Figure 4.1: Mesh for option of Table O with K = Ki = 100 



for the caU options for the two sphtting strategies. In Figure SHI we plot these meshes for the first 
example with parameters {K^ A) — {Ki, ^2)- We observe that the refinement is done mainly near 
the interface separating the region where the function vanishes and the region where it is positive. 
On these figures, we can see that the GRS is more robust and more able to follow the interface 
than the FAS. 

4.1.2. Examples in dimension 3 and 4. In this paragraph, we consider two examples of 
basket options on independent assets: one example in dimension 3 and one in dimension 4. We 
will run each example with two different strike prices. 

We give in Table the values of the truncated estimations for Ai = 12 and A2 = 13 for a 
3 dimensional basket option. The indicator error for this example is based on -E'3,3,i8,24,_r(/) and 
the number of function evaluations is given by 2000 x 1592 x 3 ~ 9.5 x 10^. 





V 


U 


C 


Cold 


iKi,A,) 


14.80805242 


2.2717704262 


1 X 10-7 


4 X 10-« 


iKi,A2) 


14.80805257 


2.2717705311 


7 X 10-** 


7 X 10- ' 


iK2,A,) 


2.927052540 


16.212009773 


6 X lG-« 


8 X 10-^ 


iK2,A2) 


2.927053375 


16.212010568 


2 X 10-« 


5 X 10-'^ 



Table 4.3: d = 3, T = 3, r = 0.05, p = 0.3, S^^'' = S*^^^ = S^^^'' = 30 
= (72 = 0-3 = 0.2, Ki = 90, K2 = 120 



We observe that the GRS is a lot more accurate than the FAS. The accuracy is about 7 
digits on the at the money option and 8 on the one out of the money. On the out of the money 
option, the FAS hardly achieves 5 digits of accuracy. 

We give in Table the values of the truncated estimations for Ai — 5 and A2 = 6. The 
indicator error for this example is based on -E4.3,i8,24,_r(/) and the number of function evaluations 
is given by 2000 x 9121 x 4 ~ 4.4 x 10"^. 

On this last example, the improvement of the algorithm is even more impressive. The accuracy 
is about 7 digits with the GRS compared to hardly 4 with the FAS. The difference between the 
two values of the criterion C for the two values of A is not due to a wrong computation of the 
integrals but simply because the integration domain was truncated too much when A = Ai. 
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V 


U 


C 


Cold 




4.22830628 


0.32667437 


1 X lO"-"^ 


3 X 10-^ 


iKi,A2) 


4.22832492 


0.32667871 


1 X lO-'^ 


4 X lO-'' 




0.16841321 


5.77905377 


7 X 10"^ 


7 X 10-'' 




0.16842047 


5.77906874 


6 X 10-*^ 


6 X lO-'' 



Table 4.4: d 4, T = 1, r = 0.05, p = 0.3, s'^^ ^ S^^'> = S*^^^ = 5^^^ = 20 
CTl = CT2 = 0-3 = cr4 = 0.1, Ki = 80, K2 = 90 



4.2. Put on minimum options. For the pricing of basket options, our approximation 
criterion was relying on the parity call put formula which does not hold for on minimum or digital 
options. It required also the pricing of a call and a put to be computed. To define error criteria 
for general options, we exploit the stochastic nature of the GRS. Instead of running only once this 
algorithm, we run it ten times only on the call pricing and compute its average V, its empirical 
variance and our error criterion will be simply Err = s^,. We also give the approximation Void 
obtained with the FAS and the estimation Vmc obtained using the crude Monte Carlo method 
for the computation of the price as the expectation of a normal random vector as described in 
Section [2m This Monte Carlo estimation is accurate up to three digits in relative error and we 
also give the number of samples required to reach such an accuracy. 

4.2.1. Examples in dimension 2. Numerical results are given in Table [53] for both cases 
Ai = 12 and A2 15. 





V 


Void 


Vmc 


Err 


Exl, Ai 


2.10306340730 


2.10306339974 


2.104291 


1.5 X 10-^" 


Exl, A2 


2.10306340508 


2.10306346643 


2.104291 


1.2 X 10-1" 


Ex2, Ai 


6.32237986596 


6.32237987060 


6.325378 


2.2 X 10-"' 


Ex2, A2 


6.32237986541 


6.32237985738 


6.325378 


1.2 X 10-1" 



Table 4.5: Put on minimum options in dimension d = 2 

Exl: T 1, r = 0.05, S^^^ = 5^ ^ = 50, p = 0.1, K ^ 45, cti = erg = 0.2 

Ex2: T = 1, r = 0.05, 5^^^ = 5^^^ = 50, p = 0.9, K = 55, cti = 0-3 = 0.2 



We observe that the GRS provides an estimation of the price which is accurate up to 10 digits. 
Indeed the integrals for the two values of A are computed with an accuracy of more than 10 digits 
and they have 10 common digits. The estimations obtained with the FAS are slightly less accurate 
up to 8 digits. In Figure [JUl we plot the two meshes for Exl in the case A — Ai. We observe 
on the GRS that the refinement is still done near the interface separating the region where the 
function vanishes and the region where it is positive but also near a line in the positive part. This 
additional line is explained by a lack of regularity of the integrand due to the minimum function 
in its definition. On these two meshes, we can see that once again the GRS is better than the FAS 
to follow the interfaces where the regularity changes. For both examples, one run of our algorithm 
requires 1.6 x 10^ function evaluations meanwhile there are 4.6 x 10^ and 1.6 x 10^ Monte Carlo 
samples for Exl and Ex2, respectively. Even if we take into account that 10 runs are used in 
order to calculate the error indicator, the number of function evaluations for both methods (GRS 
and crude Monte Carlo) are quite similar. This remains true in the following examples in higher 
dimensions. 

4.2.2. Examples in dimensions 3 and 4. We give in Table [1?^ the values of the truncated 
estimations for Ai = 12 and A2 = 15 in dimension 3. 
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(a) Geometrical Random Splitting (b) Fully Adaptive Splitting 

Figure 4.2: Mesh for option of Table [13] with if = 45 (Exl) 





V 


Void 


Vmc 


Err 


Ex3, Ai 


2.89538461 


3.021407 


2.898180 


6.3 X 10-« 


Ex3, As 


2.89538389 


2.909697 


2.898180 


3.1 X 10-« 


Ex4, Ai 


6.85473710 


6.802445 


6.854480 


6.3 X 10-« 


Ex4, A2 


6.85473692 


6.795585 


6.854480 


6.3 X 10-« 



Table 4.6: Results for put on minimum options in dimension rf = 3 
with r = 0.05, s'^'^ = S^^^ = S^^^ = 50, T = 1 and ai = a2 = crs ^ 0.2. 
Ex3: p = 0.1, if = 45 
Ex4: p = 0.9,ii: = 55 



We can see that the GRS is still very efficient and achieves an accuracy of about 7 digits. 
Now, the FAS gives poorly accurate results and is also clearly less efficient than crude Monte 
Carlo integration. This conffims the difficulties for the FAS to capture the interface when the 
dimension of the problem increases. We give in Table H771 the values of the truncated estimations 
for Ai = 12 and A2 = 15 in dimension 4. The conclusions are the same than in dimension 3. 
The FAS fails to converge while the GRS is still very accurate. On all examples, whatever the 
dimension is (up to 4), this latter strategy outperforms the crude Monte Carlo method as we have 
at least 6 digits of accuracy instead of 3 for a similar complexity. 





V 


Void 


Vmc 


Err 


Ex5, Ai 


3.567971 


3.560892 


3.574086 


6.3 X 10-^ 


Ex5, A2 


3.567971 


3.322140 


3.574086 


3.1 X 10-'' 


Ex6, Ai 


7.212993 


6.693007 


7.215822 


3.1 X 10-^ 


Ex6, A2 


7.212994 


7.305357 


7.215822 


3.1 X 10-^ 



Table 4.7: Results on minimum options in dimension rf = 4 

with r = 0.05, 5^^^ = S*^^^ = S*^^^ = 5^''' = 50, T = 1 and tji = (73 = ^3 = 

CT4 = 0.2 

Ex5: /9 = 0.1, if = 45 
Ex6: p = 0.9,is: = 55 
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4.3. Digital options. As we will see in the following numerical examples, the pricing of 
digital options is somehow a much harder problem than the two previous ones. In fact, it is close 
to the problem of looking for a small subdomain (where the function is positive) somewhere in a 
large hyperrectangle. This is well illustrated in Figure and even more in Figure [53] where this 
small subdomain is delimited by three lines. Especially during the first iterations of the algorithm, 
it may happen that a hyperrectangle of our mesh contains a part of this subdomain but with no 
quadrature points lying in it. In this case, our error indicator is zero because the function is 
considered as the null function and consequently this hyperrectangle is never split again. For 
basket or put on minimum options, we had been able to handle this sampling problem by putting 
additional quadrature points in each corner of the hyperrectangle. This trick is inefficient in the 
case of digital options because the subdomain with positive values is small and has a completely 
different form. In order to solve at least partially these convergence problems, we have increased 
the parameter a in the number of quadrature points M — ax to obtain a larger probability 

to detect a non-zero value of the function inside the hyperrectangles. Finally, we also compute 
the median Med{V) of 10 runs of the algorithm because it is a more robust estimator in case of 
false convergence. 



4.3.1. Examples in dimension 2. Our numerical results for both cases Ai = 12 and 
yl2 = 15 are given in Table 





V 


Med{V) 


Vmc 


Err 


Ex7, Ai, a = 3 


2.30072052 


2.30072041 


2.299709 


3.1 X 10-"^ 


Ex7, A2, a = 3 


2.30071826 


2.30071825 


2.299709 


3.1 X 10-^ 


Ex8, Ai,a = 3 


0.12540527 


0.15675651 


0.15600 


6.3 X 10-^ 


Ex8, A2, 


0.13848118 


0.15675549 


0.15600 


2.8 X 10-^ 


Ex8, Ai, a = 15 


0.15693827 


0.15693825 


0.15600 


3.1 X 10-^ 


ExS, A2, a = 15 


0.15681002 


0.15675531 


0.15600 


9.5 X 10"^ 



Table 4.8: Results for digital call options in dimension d = 2 
with T = 1, r = 0.05, S^^^ = S^^^ = 50, Ui = U2 = 60, cti = da = 0.2 
Ex7: p = 0.1, if = 45 
ExS: p = 0.9, if = 55 



On Ex7, we can see that we have no problem of convergence even with a — 3. The price is 
approximately equal to 2.300718, mean and median are at least 6 digits close in all cases which is 
confirmed by the error indicator. We also have 3 common digits with the Monte Carlo estimator. 
The mesh obtained in Figure 14.31 shows that the subdomain where the function is positive is 
roughly a small triangle with a surface ^ lower than 8. The ratio Ag/^ ~ 28 is still small enough 
to avoid sampling problems. 
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(a) Geometrical Random Splitting, Ai 

Figure 4.3: Mesh for option Ex7 of Table [i?^ a = 3 





-0.5 0.5 



1.5 2 2.5 



(a) Geometrical Random Splitting, A2 (b) Geometrical Random Splitting, Ai 

Figure 4.4: Mesh for option Ex8 of TableliTHla = 3 



The situation is completely different for Ex8 when a = 3. Mean and median have no more 
than 2 common digits, the error indicator is also close to 0.01. However the median has 3 common 
digits with the crude Monte Carlo estimator which means that the big convergence problems do 
not happen so often. We plot in Figure [13] two examples of meshes obtained with Ai and A2 in 
situations where the algorithm did not converge. We observe that the refinement is incomplete 
in both cases and that the one obtained with A2 misses half of the region of interest. If we add 
quadrature points by taking a — 15, we recover a high accuracy of 7 digits on the integrals at least 
when Ai = 12. The subdomain we obtain in this case in Figure [53] has now a triangular shape. 
Its surface is now roughly equal to 1 which explains the convergence problems we encountered 
when a = 3 as the ratio A^/S, ~ 225. 
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-3 

-3-2-101234 

(a) Geometrical Random Splitting, Ai 

Figure 4.5: Mesh for option Ex8 of Table [121 « = 15 



For a = 15, one run of the GRS requires 8 x 10^ function evaluations, which is 5 times more 
than in the case a — 3 meanwhile there are 4.4 x lO"^ and 1.4 x 10^ Monte Carlo samples for 
respectively Ex7 and ExS. 

4.3.2. Examples in dimension 3. We give in Table [121 the values of the truncated esti- 
mations for A2 — 12 and A2 = 15. 





V 


Med{V) 


Vmc 


Err 


Ex9, Ai, a = 20 


1.64950182 


1.64950187 


1.646800 


2.5 X 10"*^ 


Ex9, A2, a = 20 


1.64948232 


1.64948236 


1.646800 


1.9 X 10"*^ 


ExlO, Ai, a = 40 


0.09316076 


0.09316072 


0.093269 


3.1 X 10"^ 


ExlO, ^2, a = 40 


0.09307638 


0.09316133 


0.093269 


2.5 X 10-4 



Table 4.9: Results on digital options in dimension d = 3 

with r = 0.05, T = 1, 5^^^ = S^^^ = S^^^ = 50, Ui = U2 = U3 = 60 

Ex9: p = 0.1,K = 4:5 

ExlO: p = 0.9,X = 55 



The number of function evaluations are 6.3 x 10^ for a = 20 and 1.3 x 10* for a = 40 meanwhile 
there are 3.0 x 10^ and 7.3 x 10^ Monte Carlo runs respectively for Ex9 and ExlO. On Ex9, we 
need to increase a to converge up to 5 or 6 digits but the number of function evaluations is still 
comparable to the size of the Monte Carlo sample. On ExlO with A2, the Monte Carlo method has 
a smaller variance and the GRS a weaker accuracy which makes these two methods comparable in 
terms of efficiency. However, with Ai, the results obtained with the GRS seems very accurate in 
comparison with those obtained with the Monte Carlo Method. A more detailed discussion about 
the influence of the initial volume of the hypperrectangle will be given in dimension 4. 

4.3.3. Examples in dimension 4. We give in Table 14.101 the values of the truncated 
estimations for the usual values Ai — 12, A2 — 15 and also for two smaller values ^3 = 5 and 
A4 = 6. 
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Med{V) 


Vmc 


Err 


Exll, Ai, a = 30 


1.22934667 


1.22934613 


1.228935 


3.1 X 10"*^ 


Exll, A2, a = 30 


1.22934272 


1.22934341 


1.228935 


3.1 X 10"" 


Exl2, Ai, a = 40 


0.04827239 


0.04826375 


0.060981 


1.2 X 10"^ 


Exl2, A2, a 40 


0.04056702 


0.03601258 


0.060981 


9.5 X 10"-^ 


Exl2, A3, a = 40 


0.06126412 


0.06126407 


0.060981 


6.3 X 10-^ 


Exl2, yl4, a = 40 


0.06126593 


0.06126569 


0.060981 


6.3 X 10- ' 



Tabic 4.10: Results on digital options for d = 4 

with r = 0.05, T=l, S^a^ = S^^^ = 5^^^ = 5^"*^ = 50, Ui =^ U2 ^ U3 = Ui ^ 60 

Exll: p = 0.1,ii: = 45 

Exl2: r = 0.05, p = 0.9, if = 55 



The number of function evaluations are 4.4 x 10^ for a — 30 and 5.8 x 10* for a = 40 meanwhile 
there are 2.2 x 10^ and 4.3 x 10" Monte Carlo runs respectively for Exll and Exl2. Obviously, 
the complexity of the GRS increases with the dimension but it is still more efficient on Exll than 
the Monte Carlo method. On Exl2, the volume of the subdomain where the function is positive is 
getting too small compared to the size A| of the initial hypperrectangle for the GRS to converge. 
Nevertheless, we have been able to recover a very good accuracy on Exl2 by drastically reducing 
the size of the initial hypercube using A3 and A4. Indeed, the corresponding price values have 
almost 6 common digits and the error criterion is 2 x 10~^. This shows that is would be worth 
using some basic preliminary numerical methods to circumvent the region of interest in a relatively 
small box. 

5. Delta computation. 

5.1. Tchebychef Interpolation and differentiation. To compute the partial derivatives 
involved in the computation of the Delta, we use a standard approach based on Tchebychef in- 
terpolation polynomials. We describe in detail the two dimensional case for a general function 
/(x, y). First, we compute the Tchebychef interpolation polynomial Tm,h,XQ.yo i^) of /(x, yo) using 
m interpolation points in the interval ]x — x + h[. The approximation of the first component 

d 

of the gradient of / at point {xo,yo) is given by -g-T,n,h,xo,yo{^)- Thus, it requires to compute 
m prices, which is done using the GRS. Using different values of the parameters h and m, we are 
able to obtain reliable Delta approximations. 

5.2. Monte Carlo approach. To measure the efficiency of the approach described in Para- 
graph 15.11 we need an alternative method to compute the Delta, which is sufficiently accurate 
to serve as a benchmark. We relied on a finite difference approach coupled with a Monte Carlo 
technique. We closely followed the recommendations of [5] to tune the number of Monte Carlo 
simulation and the finite difference step accordingly to ensure the best possible accuracy. For a 
Monte Carlo method with n samples, we used a finite difference step equal to /i„ = n~^^^ so that 
under mild assumptions the finite difference estimator converges almost surely to the exact value 
and satisfies a central limit theorem with the rate n"^/'^. 

5.3. Numerical Results. For symmetry reasons, we compute only the first component of 
the Delta. We compare in table l5. li the Monte Carlo approach and the method based on Tchebychef 
interpolation on 4 examples. These latter are taken in dimension 3 and 4 for basket and minimum 
options. For the interpolation method, we use two sets of parameters chosen to have two different 
accurate estimations of the Delta. The maximum number of function evaluations required for the 
examples in dimension 4 is 5 x 4.4 x 10^ = 2.2 x 10* which is comparable to the number of samples 
n used in the Monte Carlo method. 
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Monte Carlo 


m = 3,h = 0.05 


m = 5,h = 0.1 


n 


Exl3 


0.300088 


0.3002853 


0.3002864 


1.4 X 10** 


Exl4 


-0.240865 


-0.2382143 


-0.2382098 


1.6 X 10« 


Exl5 


0.230224 


0.2303219 


0.2303214 


1.1 X 10** 


Exl6 


-0.186628 


-0.1837101 


-0.1836998 


1.6 X 10** 



Table 5.1: Delta Computations with T=l,r = 0.05, a = 0.2, = 50 

Exl3: Basket options, d = 3,p = 0.1, K ^ 45 

Exl4: Put on minimum options, d = 3, p = 0.5, K = 55 

Exl5: Basket options, d = 4,p = 0.1, K = 45 

Exl6: Put on minimum options, d = A,p = 0.5, K = 55 



The value of n is chosen to ensure an accuracy of 3 digits on the Delta. All estimations 
of the Delta using the interpolation method have about 3 common digits with the Monte Carlo 
estimation. Since for the two sets of parameters we have obtained at least 5 common digits for the 
two different estimations, we may conclude that the interpolation method is more efficient than 

the Monte Carlo approach. 

6. Dimension reduction and control variates for high dimensional problems. 

6.1. Description of the method. The GRS developed and tested in the previous sections 
suffers from the curse of dimensionality but shows a very impressive accuracy in low or medium 

dimensions up to 4 or 5. In this section, we propose a way to use this method to devise a control 
variate for high dimensional models. We are still interested in computing expectations of the form 
E(^(5t)) where St is defined as in the previous section. The basic idea of the method is to 
perform a principal component analysis of St in order to reduce the dimension of the problem 
by keeping only the leading components and setting the others to zero. The expectation in the 
reduced model can be computed quickly and accurately using the GRS and can serve as a control 
variate for the original problem. 

6.1.1. Principal Component Analysis. We rewrite the model to embed the correlation 
matrix F into the volatility structure turned into a matrix S defined by 

S = diag(c7i, ...,(7d) r diag((Ji, ...,(Td). 

Note that the vector (ctiWt, ...,(J2Wt) is a Gaussian vector with covariance matrix a/TS. It is 
straightforward to check that the matrix S inherits its symmetric positive definite feature from 
the one of T and hence admits an orthonormal basis of eigenvectors with positive eigenvalues. 
Let D be the diagonal matrix built up with these eigenvalues sorted in decreasing order and 
let P be the corresponding matrix of eigenvectors sorted accordingly. Note that reordering the 
columns of P does not change its orthonormal property and that wc have P~^ = P* and as a 
consequence S = P*DP. Now, if we let = diag(V-Dii, \/Dm) and define the symmetric 
matrix H = P^D^P, we obtain S = H*H. This leads to the identity in distribution 

{giWt,...,<JiWt) = VtHG 
where G is a standard normal vector in M**. Then, we can write the identity in law 

S't = Si, exp (^{r - + VfHiG^ 

where Hi denotes the i — th row of H. This new expression for S can help us to devise a reduced 
model S with effective dimension I < d such that S can be written 

4 = Si exp (r-^)T+ y/fHid\ 



withG = (Gi,...,G(,0,...,0). 
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6.1.2. Control variates. If I is sufficiently small, namely less than 3 or 4, we can compute an 
approximation I of E ('(/;( St)) using the GRS both very accurately and with a small computational 
cost. If the dimension reduction works well, we can also expect the random variable V'(*S't) to be 
close to iplSr) and as a consequence to appear as a natural choice for a control variate. Our new 
Monte Carlo estimator for K{i/j{St)) will be 

-, Nmc 

(1) (i) 

where the Nmc random samples Srf and Sj! are built using the same independent drawings 
G(J) from the standard normal distribution over M . Because the effective dimension Z of G is a 
lot smaller than d and the GRS is so efficient in low dimensions, one can implement this control 
variate approach with almost no extra computational cost. The idea of using such a control variate 
was already proposed for example in 4 but it was limited to a sum of control variates depending 
on only one variable. Indeed this computation was relying on closed forms for the expectations 
of the control variates which are only available in dimension one. Hence, our new approach is far 
more general. 

6.2. Numerical results. 

6.2.1. Toy correlation structures. All the examples deal with basket options in dimension 
5 with three different values of the correlation coefficient p from very correlated assets {p = 0.9) to 
negatively correlated assets (p = —0.1). The different estimations of the price and its confidence 
interval, named Ieq , Iei , ■, , using the control variate method are given in table 16.11 We 
compute the estimation of the price /s using the GRS in order to obtain a reference value. 





Ieo 


Iei 


Ie2 


lE:-i 


h 


Exl7 


8.61236 ±0.020 


8.61333 ±0.0013 


8.61357 ±0.0010 


8.61407 ±0.0003 


8.61404 


Exl8 


7.51683 ±0.0130 


7.52217 ±0.0072 


7.52395 ±0.0042 


7.52621 ±0.0023 


7.52490 


Exl9 


7.29012 ±0.0103 


7.28436 ± 0.0074 


7.27969 ±0.0042 


7.27931 ±0.0038 


7.27548 



Table 6.1: Basket options with 0.95% confidence interval 

= 5, T = 1, r = 0.05, K = 45, So = 50 
a = (0.156,0.442,0.325,0.134,0.114) and 100,000 samples 
Exl7: p = 0.9, Exl8: p = 0.1, Exl9: p = -0.1 



We observe that the control variate technique gives very good results. Especially with 3 PCA 
components, the results are very close to I^. Moreover, the GRS becoming very costly in dimension 
5, the control variate approach with 3 components seems now more efficient. 

As it might be expected, we also noticed that when the assets are very correlated the variance 
reduction is more efficient. With p — 0.9, the first component of the PCA contains almost all the 
information of the model which ensures a huge variance reduction. When p = —0.1, even with 
3 PCA components, the variance reduction is so small that it hardly compensates the additional 
cost of the control. Indeed this latter is twice bigger than the initial Monte Carlo cost. 

6.2.2. A more general correlation structure. In this paragraph, we are concerned with 
testing our approach on more realistic covariance structures. We consider a block covariance 
matrix Gex20, in which the assets can be split into two subsets which are negatively correlated, 
whereas inside each set all the assets are positively correlated with the same correlation. Such 
covariance structures are quite common in practice. We can find in Table \^7I\ the results obtained 
with a standard Monte Carlo and our control variate approach method for different numbers of 
PCA components. At first sight, one could think that such a covariance structure could cause 
trouble to our approach and yet it manages to divide the width of the confidence interval by a 
factor of 10, which is really impressive if one remembers that to achieve such a confidence interval 
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with a crude Monte Carlo it would require 100 times more samples. We can see on this example 

that to ofFcctivcly reduce the variance, wc need at least two PCA components and no closed form 
formulae arc available for basket options in dimension 2 and 3. This example highlights the key 
role played by the GRS approach. 



Ex20 



( 1 
0.8 
0.8 
0.8 
0.8 
-0.5 
-0.5 
-0.5 
-0.5 
\ -0.5 



0.8 
1 

0.8 
0.8 
0.8 
-0.5 
-0.5 
-0.5 
-0.5 
-0.5 



0.8 
0.8 
1 

0.8 
0.8 
-0.5 
-0.5 
-0.5 
-0.5 
-0.5 



0.8 
0.8 
0.8 

1 

0.8 
-0.5 
-0.5 
-0.5 
-0.5 
-0.5 



0.8 
0.8 
0.8 
0.8 
1 

-0.5 
-0.5 
-0.5 
-0.5 
-0.5 



-0.5 
-0.5 
-0.5 
-0.5 
-0.5 
1 

0.4 
0.4 
0.4 
0.4 



-0.5 
-0.5 
-0.5 
-0.5 
-0.5 
0.4 
1 
0.4 
0.4 
0.4 



-0.5 
-0.5 
-0.5 
-0.5 
-0.5 
0.4 
0.4 
1 

0.4 
0.4 



-0.5 
-0.5 
-0.5 
-0.5 
-0.5 
0.4 
0.4 
0.4 
1 

0.4 



-0.5 
-0.5 
-0.5 
-0.5 
-0.5 
0.4 
0.4 
0.4 
0.4 
1 







Si 


£^2 


^3 


Ex20 


3.1912 ±0.011 


3.1899 ±0.009 


3.1908 ±0.002 


3.1906 ±0.001 



Table 6.2: Prices with 0.95% confidence interval for homogeneous call Basket 
options on example Ex20 with 

T = 2,r = 0.02, K = 105, a = 0.2, = 100 and 100, 000 samples 



7. Conclusion. In this paper, we have developed a new numerical integration algorithm 

based on a Geometrical Random Splitting. This new algorithm has been successfully applied to 
the pricing of Vanilla options in dimensions up to five. In particular, this new algorithm efficiently 
handles the difficult problems of pricing and hedging digital options. In most situations, the 
accuracies we have obtained were out of reach for a crude Monte Carlo approach. For higher 
dimensions, we have shown that the GRS can be used as a control variate when associated to a 
principal component analysis. The resulting variance reduction was quite impressive especially 
for very correlated assets. One important issue would be the development of other dimension 
reduction techniques that could be coupled with the GRS. The GRS itself can obviously have 
many other applications for instance in computer experiment or in the adaptive approximation of 
partial differential equations. The GRS algorithm proposed in this paper is very satisfactory on 
most examples studied. 
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